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1 Introduction 

Irrigation  pumping  is  concentrated  during  a three-  to  five-month  period  for  much  of  the  pump- 
irrigated  areas  of  the  U.S.  Electric  powered  pumping  increases  the  peak  demand  of  the  electric 
utilities  serving  these  irrigation  loads  whether  they  are  wholesale  or  retail  loads.  Although  pumping 
increases  peak  demands  for  the  electric  power  supply  system,  it  is  also  a significant  part  of  the  annual 
revenue  for  some  utilities  and  contributes  to  the  economies  of  communities  in  irrigated  areas. 

1.1  Objective 

Our  objective  was  to  evaluate  or  develop  electrical  load  prediction  programs  and  irrigation  schedul- 
ing programs  that  predict  irrigation  requirements  for  integration  into  an  interactive  load  man- 
agement program  that  could  be  used  by  electric  utilities  for  controlling  peak  demands  where  an 
irrigation  load  is  a significant  part  of  the  demand. 

1.2  Procedure 

The  procedure  was  to  review  the  literature,  evaluate  the  AGEND  Irrigation  Model  for  forecasting 
irrigation  electrical  use,  establish  and  meet  with  a customer  focus  group  and  acquire  data  to  begin 
development  of  a forecasting  model. 

Many  of  the  attached  references  were  well  known  to  the  principal  investigators.  They  were 
used  for  the  AGEND  model  and  several  studies  completed  under  the  guidance  of  the  principal 
investigator.  The  later  studies  were  conducted  in  the  1970’s  and  studied  the  load  management 
from  the  irrigator’s  viewpoint.  The  current  study  was  focused  at  the  level  of  the  supplier  to 
the  irrigator.  The  forecasting  model  that  was  evaluated  and  developed  will  compliment  previous 
studies.  A number  of  references  identified  in  the  automated  literature  search  were  written  by  the 
principal  investigators  and  are  cited  here. 

The  AGEND  Irrigation  Model  was  developed  for  EPRI  under  a contract  with  Quantum  Con- 
sulting, Inc.  AGEND  is  a detailed  model,  but  when  given  the  necessary  input  data  it  should  provide 
a good  forecast  of  electrical  demand.  The  model  has  not  been  readily  accepted  by  the  utilities. 
The  discussion  with  potential  users  with  the  customer  focus  group  showed  that  the  model  required 
more  detailed  data  than  is  readily  available.  The  model  required  that  the  crop  must  be  known  for 
each  field  irrigated.  In  addition,  the  data  required  for  each  system  included  pump  efficiency,  motor 
size,  pumping  lift,  pipe  sizes  and  discharge.  The  estimated  running  time  is  generated  from  crop 
water  requirements.  Though  the  model  can  give  accurate  electrical  usage,  obtaining  the  detailed 
information  necessary  to  run  the  model  is  difficult.  > 

A customer  focus  group  made  up  of  TRI-STATE  G&T  technical  representatives  and  most  of  the 
REA  rural  managers  and  technical  specialists  with  significant  irrigation  loads  in  eastern  Colorado 
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and  western  Nebraska  provided  input  on  their  needs  for  forecasting  electrical  loads.  They  indicated 
the  type  of  data  that  are  available  for  use  in  a model.  Several  cooperatives  are  either  considering 
or  are  starting  load  control  programs  for  irrigation  loads.  They  stated  that  their  need  is  for  short 
term  forecasting  of  irrigation  loads  with  emphasis  on  the  shoulder  months  of  June  and  September. 
Just  forecasting  the  next  couple  of  days  would  be  valuable  for  use  in  load  control  programs.  This 
would  allow  control  to  be  minimized  and  provide  the  most  benefit  for  their  program. 

They  identified  the  types  of  data  that  are  readily  available  through  either  Coop  offices  or 
from  TRI-STATE.  The  available  data  varied  from  Coop  to  Coop  but  usually  included  billing  point 
demand  and  total  kWh  delivered.  Sometimes  the  crop,  demand,  kWh,  location,  water  use  and  type 
of  irrigation  system  are  available.  The  customers  were  all  willing  to  share  data  for  our  analysis.  We 
have  obtained  the  demand  data  from  all  billiing  points  from  Highline  Electric  for  our  first  analysis. 
The  service  areas  also  have  weather  station  data  available  which  will  allow  calculation  of  crop  water 
use. 

Our  initial  approach  was  to  perform  a time  series  analysis  of  the  demand  data  for  the  years 
1992-1996.  The  S-Plus  statistical  package  was  purchased  for  doing  this  part  of  the  analysis.  The 
objective  was  to  first  determine  the  maximum  information  that  could  be  obtained  from  the  data 
that  was  readily  available  and  was  already  computer  compatible.  The  data  were  obtained  directly 
from  TRI-STATE  which  had  edited  the  data  for  quality  assurance. 

We  outlined  our  approach  to  this  study  to  TRI-STATE  personnel  and  obtained  their  ideas  as 
to  the  needs  and  benefits  for  load  management.  We  need  to  understand  not  only  the  needs  of  the 
rural  cooperatives  and  their  constraints  but  also  the  objectives  and  approach  of  the  supplier  to 
deliver  the  necessary  power  to  meet  demand.  It  was  recognized  that  though  each  has  the  objective 
of  reducing  costs,  the  paths  to  solution  can  often  lead  to  conflict.  It  confirms  the  benefits  for  the 
rural  cooperatives  for  a better  procedure  to  forecast  short  term  demands.  The  electrical  delivery 
system  is  one  where  each  entity  is  optimizing  within  the  constraints  of  their  own  system  and  the 
operating  procedures  of  the  other. 

The  demand  data  from  Highline  Electric  were  obtained  for  the  years  1992-1996.  Climatic  data 
from  the  area  were  used  to  estimate  crop  water  use  and  rainfall  amounts.  The  initial  analysis 
was  to  develop  correlations  between  the  daily  peak  demands  and  the  climatic  data.  The  primary 
correlation  was  determined  between  the  rainfall  and  sharp  drop  in  peak  demands  from  one  day  to 
the  next.  As  interesting  as  this  may  be,  it  did  not  provide  the  needed  estimate  of  the  future  peak 
daily  demands.  The  relationship  between  peak  demands  and  ETA  was  not  found  to  be  significant. 
It  may  be  that  the  irrigators  did  not  schedule  irrigations  based  on  crop  ET  but  based  their  decisions 
on  perception  and  historical  experience  for  scheduling  irrigations. 

The  remainder  of  this  report  discusses  our  efforts  to  create  an  algorithm  for  predicting  the 
maximum  daily  power  usage  which  could  be  implemented  with  ease.  It  was  desired  that  this 
algorithm  be  capable  of  making  predictions  several  days  in  advance.  In  particular,  seven-step 
predictions  were  of  interest  in  order  to  allow  the  involved  parties  sufficient  time  to  react. 

The  data  analyzed  in  detail  was  the  power  demand  data  (in  kilowatts)  recorded  in  thirty  minute 
intervals  at  the  Holyoke  station  for  the  months  of  April  through  October,  years  1992  to  1996. 
Additionally,  daily  weather  data  in  the  general  vicinity  of  the  Holyoke  station  for  the  months  of 
April  through  September  were  included  in  the  analysis.  The  climatic  data  consisted  of  maximum 
daily  temperature,  minimum  daily  temperature,  daily  precipitation  (in  mm),  ETP  (as  calculated 
by  the  Penman  equation),  and  ETAcorn  (a  variation  of  ETP  adjusted  for  corn,  the  major  crop  in 
the  area). 
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2 Exploratory  Plots 

2.1  Condensing  the  Original  Power  Usage  Data 

The  Holyoke  power  data  were  first  condensed  by  finding  the  maximum  power  usage  for  each  day 
and  discarding  the  remaining  points.  Since  our  goal  was  to  predict  maximum  daily  power  usage, 
no  information  seemed  to  be  lost  in  using  these  subsets  of  daily  maxima  as  our  working  data  sets. 
This  was  apparent  when  comparing  time  series  plots  of  the  original  data  and  the  daily  maximum 
data  subset  depicted  in  Figures  1 through  10.  Additionally,  we  only  considered  the  power  data 
through  the  month  of  September  since  most  irrigations  are  completed  by  mid  September. 

2.2  Data  Exploration 

We  examined  several  time  series  plots  and  scatter  plots  in  order  to  obtain  a better  understanding 
of  how  peak  usage  changed  over  time  and  what  were  some  of  the  factors  influencing  this  movement. 
Initially,  we  considered  the  relative  differenced  data  Xt  (Xt  = ^yJYt , where  Yt  is  the  maximum 
daily  power  usage)  in  an  effort  to  create  a more  stationary  sequence.  This  series  is  depicted  for 
all  available  years  in  Figures  11  through  15.  The  darkened  circles  indicate  that  the  precipitation 
on  day  t was  positive  (Rt  > 0 ram).  In  studying  these  graphs,  it  appears  that  there  is  an  inverse 
relationship  between  rainfall  and  power  usage.  This  relationship  is  depicted  more  clearly  in  the 
scatter  plots  of  Xt  vs.  Rt  (Figures  16  through  20).  A regression  line  was  added  to  these  plots  for 
clarification.  Although  the  strength  of  the  relationship  between  these  two  variable  changes  from 
one  year  to  the  next,  it  is  evident  in  all  five  years. 

Scatter  plots  involving  the  remaining  weather  variables  did  not  reveal  any  further  relationships. 
In  addition  to  our  consideration  of  various  scatter  plots,  we  also  examined  several  brush  plots  in  the 
hopes  of  uncovering  further  associations  among  the  variables.  A brush  plot  is  an  S-Plus  graphing 
device  that  allows  the  user  to  interact  with  various  2-dimensional  and  3-dimensional  scatter  plots. 
One  of  the  options  allows  the  user  to  highlight  specific  points  on  one  graph  resulting  in  these  points 
being  highlighted  on  all  of  the  remaining  graphs.  This  is  often  useful  for  discovering  relationships 
which  are  only  present  for  a subset  of  the  data.  The  user  may  also  change  the  variable  axes 
on  the  3-dimensional  graph  to  explore  relationships  between  any  three  variables  and  rotate  this 
graph  to  obtain  different  perspectives.  However,  no  further  relationships  were  discovered  using  this 
procedure.  Our  next  step  was  to  begin  modeling  the  series. 
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Figure  1:  Time  Series  Plot  of  Power  Usage  Data  for  1992 


Figure  2:  Time  Series  Plot  of  Maximum  Daily  Power  Usage  for  1992 
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Figure  3:  Time  Series  Plot  of  Power  Usage  Data  for  1993 


Figure  4:  Time  Series  Plot  of  Maximum  Daily  Power  Usage  for  1993 
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Figure  5:  Time  Series  Plot  of  Power  Usage  Data  for  1994 
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Figure  6:  Time  Series  Plot  of  Maximum  Daily  Power  Usage  for  1994 
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Figure  7:  Time  Series  Plot  of  Power  Usage  Data  for  1995 


Figure  8:  Time  Series  Plot  of  Maximum  Daily  Power  Usage  for  1995 
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Figure  9:  Time  Series  Plot  of  Power  Usage  Data  for  1996 
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Figure  10:  Time  Series  Plot  of  Maximum  Daily  Power  Usage  for  1996 
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Figure  11:  Time  Series  Plot  of  Relative  Differenced  Data  for  1992 


Figure  12:  Time  Series  Plot  of  Relative  Differenced  Data  for  1993 
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Figure  13:  Time  Series  Plot  of  Relative  Differenced  Data  for  1994 


Figure  14:  Time  Series  Plot  of  Relative  Differenced  Data  for  1995 
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Figure  15:  Time  Series  Plot  of  Relative  Differenced  Data  for  1996 


Figure  16:  Scatter  Plot  of  Xt  versus  Rt  for  1992 
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Figure  17:  Scatter  Plot  of  Xt  versus  Rt  for  1993 


Figure  18:  Scatter  Plot  of  Xt  versus  Rt  for  1994 
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Figure  19:  Scatter  Plot  of  Xt  versus  Rt  for  1995 


Figure  20:  Scatter  Plot  of  Xt  versus  Rt  for  1996 
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3 Initial  Model  Considerations 

3.1  Structural  Models 


In  order  to  model  the  relative  differenced  series  we  began  with  a structural  time  series  model  written 
in  state  space  form.  This  type  of  model  is  often  thought  to  be  a good  starting  point  because  it 
encompasses  a wide  range  of  time  series  models,  and  with  the  aid  of  the  Kalman  recursions,  the 
analysis  is  relatively  straightforward.  The  general  form  of  a state  space  model  for  a iu-dimensional 
time  series  {Y(,t  = 1, 2, . . .}  is  given  below. 

Yt  = GtXt  + Wt,  W*  ~ WN(0,  {Ft}) 

X<+i  = FtXt  + Vt)  Vt~WN(0,{Q<}) 

where  X*  is  a u-dimensional  state  variable,  {Gt}  is  a sequence  of  wxv  matrices,  {Ft}  is  a sequence 
of  vxv  matrices,  {V*}  is  uncorrelated  with  {Wt}  and  Xi  is  uncorrelated  with  both  {Vt}  and 
{W<}.  The  first  equation  is  called  the  observation  equation  because  it  expresses  the  observation 
vector,  Yt,  as  a linear  function  of  the  state  vector,  X*.  The  second  equation  is  known  as  the  state 
equation  because  it  defines  the  state  vector  at  time  t -f  1 in  terms  of  the  previous  state  vector  at 
time  t along  with  a noise  term. 

The  Kalman  equations  are  a recursive  set  of  formulas  which  can  be  used  to  find  the  best  linear 
mean-squared  predictor  of  the  state  vector,  X*,  based  on  the  information  available  at  time  t — 1. 
This  predictor  is  denoted  by  X*  or  Pt_x(Xt)  and  is  determined  by  the  following  equations: 

Xx  = P(Xx|Y0) 

Qx  = F[(Xx-Xi)(Xi-Xi)'] 

X(+1  = FtXt  + QtA^iYt-GtXt) 
f2t+i  = FtCltFf  +Qt  — ©(AJ-1©:, 

where 

t = 1,2,..., 

Q t = the  error  covariance  matrix  of  X< 

A t = GfQtGj  + Ft,  A^1  is  any  generalized  inverse  of  At 

©t  = fag;. 


It  also  follows  from  these  equations  that  the  best  linear  mean-squared  predictors  of  Xt+h  and  Y t+h. 
based  on  the  observations  up  to  time  t are 


FtXt+h  = [Ft+h-i.Ft+h-2  • * ’Ft+xJFtXt+X) 
FtYt+h  = Gt+hPtXt+h‘ 


The  reader  is  referred  to  Brockwell  and  Davis  [2]  for  further  information  on  structural  modeling 
and  prediction  via  the  Kalman  recursions. 

We  initially  considered  the  relative  differenced  series  to  be  a randomly  varying  trend  with 
seasonal  and  noise  components.  This  model  was  chosen  because  we  wanted  to  determine  whether 
there  was  a trend  in  the  data  and  also  whether  there  was  a weekly  effect.  The  state  space  form  for 
this  model  is  as  above  with  Yt,  Xf,  Gf,  Ft,  Rt  and  Qt  defined  as  follows: 


Yt 

Xt 


Ut+i  - Ut 


Ut  = maximum  daily  kw  usage  at  time  £, 


(Mf,  Pf , St,  St_x,  St- 2)  St_3,  St_4,  St— 5)', 
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where  M*  is  the  deterministic  signal  at  time  t,  Bt  is  the  slope  of  the  linear  trend  at  time  t and  st 
is  a seasonal  component  satisfying  st+7  = st  and  £j=1  st  = 0, 
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t = 1,2,...,  183  (time  t — 1 representing  April  1 and  time  t = 183  representing  September  30). 

This  model  was  fit  using  a Fortran  program  which  incorporates  the  use  of  Hooke  and  Jeeve’s 
algorithm  to  find  maximum  likelihood  estimators  of  o\,  and  Q via  the  Kalman  recursions.  Xi  = 0 
and  = 0 ■ I7.7  were  the  initial  conditions  used  in  the  fit.  For  further  reference  on  obtaining  these 
maximum  likelihood  estimates,  see  Brockwell  and  Davis  [2]. 

From  this  initial  fit,  the  series  did  not  appear  to  have  either  a linear  trend  or  a weekly  component 
(both  a\  and  <t|  were  approximately  0).  Therefore,  we  refit  the  model  as  a random  walk  with  a noise 
term.  This  model  fit  the  data  quite  well  considering  the  fact  that  we  had  not  yet  incorporated  any 
of  the  weather  data.  The  observed  values  (denoted  by  ‘o’)  and  one-step  predicted  values  (denoted 
by  ‘*’)  of  maximum  kw  usage  under  this  model  for  1992  through  1996  are  depicted  in  Figures  21- 
25.  The  maximum  likelihood  estimates  of  a\  and  a \ obtained  by  the  Fortran  program  are  given  in 
Table  3.1. 


Table  3.1:  Maximum  Likelihood  Estimates 


Year 

1992 

0.05 

0.07 

1993 

0.05 

0.06 

1994 

0.05 

0.03 

1995 

0.05 

0.03 

1996 

0.05 

0.06 
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o:  Y(t)  observed  *:  Y(t)  predicted  o:  Y(t)  observed  *:  Y(t)  predicted 


Figure  21:  Structural  Model  Fit  for  1992 


Figure  22:  Structural  Model  Fit  for  1993 
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o:  Y(t)  observed  *:  Y(t)  predicted  o:  Y(t)  observed  *:  Y(t)  predicted 


Figure  23:  Structural  Model  Fit  for  1994 


Figure  24:  Structural  Model  Fit  for  1995 
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Figure  25:  Structural  Model  Fit  for  1996 
3.2  Regression  Models 

In  an  effort  to  improve  upon  the  structural  model  discussed  in  Section  3.1,  we  decided  to  consider 
various  regression  models  to  determine  which  weather  variables  might  improve  our  predictions. 
One  of  the  first  regression  models  we  considered  was  the  following: 


Xt  — Po  4-  Pi  Rt  4-  j 


where: 


Xt  = 


Yt+1  - Yt 
Yt 


, Yt  — maximum  daily  kw  usage  at  time  t 


Rt  = daily  rainfall  at  time  t 
et  ~ iV(0,<72) 


We  considered  this  model  because  of  the  relationship  between  the  relative  differenced  data  and 
rainfall  evident  in  Figures  16-20.  The  ANOVA  table  values  obtained  using  S-Plus  are  given  in 
Table  3.2a. 
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Table  3.2a:  ANOVA  Table  Values 


Year 

Parameter 

Estimate 

Standard  Error 

P-value 

1992 

Po 

0.0702 

0.0214 

0.0013 

Pi 

-0.0140 

0.0029 

0.0000 

1993 

Po 

0.0850 

0.0193 

0.0000 

Pi 

-0.0217 

0.0033 

0.0000 

1994 

Po 

0.0410 

0.0159 

0.0109 

Pi 

-0.0153 

0.0041 

0.0003 

1995 

Po 

0.0234 

0.0150 

0.1194 

Pi 

-0.0058 

0.0050 

0.2422 

1996 

Po 

0.0735 

0.0225 

0.0013 

Pi 

-0.0100 

0.0024 

0.0001 

As  we  saw  in  Figures  16-20,  the  inverse  relationship  between  the  relative  differenced  data  and 
rainfall  is  quite  strong  for  all  of  the  years  except  1995. 

Using  this  regression  model  we  were  able  to  obtain  more  accurate  predictions,  (in  terms  of 
mean-squared  error)  for  the  original  series  Yt  than  we  were  under  the  structural  model  framework. 
The  root  mean-squared  errors  of  the  predictions  are  given  in  Table  3.2b  for  both  of  these  mod- 
els as  well  as  an  additional  regression  model  which  will  be  discussed  below.  The  mean-squared 
error  for  the  original  series  was  calculated  using  the  formula  MSE  = ^ 52»=i  0^t+i  ~ Yt+ 1)2  = 
y§2  SCi=i  (Yt+ 1 — Y*(At  + l))2-  The  sum  consists  of  182  terms  instead  of  183  terms  since  the  models 
were  fit  to  the  relative  differenced  series  {A*}.  The  predictions  for  the  series  {V*}  were  then  found 

J y 1 J 

from  the  predictions  of  {A*}  using  the  relationship  A*  = The  gain  in  predictive  power 

can  also  be  seen  by  comparing  Figures  21-25  with  Figures  26-30.  In  these  figures,  the  observed 
values  are  plotted  as  ‘o’  while  the  one-step  predictions  are  plotted  as  The  darkened  circles 
in  Figures  26-30  indicate  that  there  was  a positive  amount  of  precipitation  on  the  previous  day. 
Again,  the  inverse  relationship  between  power  usage  and  rainfall  is  evident. 


Table  3.2b:  Root  Mean-squared  Errors 


Model 

1992 

1993 

1994 

1995 

1996 

Model  1 

983 

1,016 

870 

745 

1,216 

Model  2 

671 

688 

810 

700 

1,152 

Model  3 

717 

753 

819 

720 

1,195 

Model  1 is  the  Structural  Model,  Model  2 is 
the  Regression  Model  with  Precipitation, 
and  Model  3 is  the  Regression  Model  with 

Precipitation  and  ETA. 
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o:Y(t)  *:  Yhat(t)  o:Y(t)  *:  Yhat(t) 

2000  4000  6000  8000  10000  12000  0 2000  4000  6000  8000  10000  12000 


Figure  26:  Regression  Model  with  Precipitation  - 1992 


Figure  27:  Regression  Model  with  Precipitation  - 1993 
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o:  Y(t)  *:  Yhat(t)  o:  Y(t)  *:  Yhat(t) 


Figure  28:  Regression  Model  with  Precipitation  - 1994 


Figure  29:  Regression  Model  with  Precipitation  - 1995 
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Figure  30:  Regression  Model  with  Precipitation  - 1996 


To  determine  whether  other  weather  variables  might  assist  in  predicting  power  usage,  we  in- 
cluded various  combinations  of  the  available  variables  as  predictors  in  a regression  model  setting. 
One  such  model  is  given  below  and  is  included  in  Table  3.2b. 


where: 


Xt  = Po  + Pi  Rt  + P2^t  + €t 


Xt 

Rt 

Et 


Yt+ 1 - Vt 
Yt 


Yt  = maximum  daily  kw  usage  at  time  t 


daily  rainfall  at  time  t 

daily  value  of  ETAcorn  at  time  t 

AffO.cr2) 


The  predictions  under  this  model  were  better  than  those  obtained  from  the  structural  model. 
However,  they  were  not  as  accurate  as  those  obtained  with  precipitation  as  the  sole  predictor 
variable  (see  Table  3.2b).  This  was  the  case  with  all  of  the  regression  models  we  considered  involving 
combinations  of  rainfall  and  other  weather  variables  as  the  explanatory  variables. 


4 Base  Model  Approach 

At  this  point,  we  began  to  re-think  our  modeling  strategy.  In  order  to  make  predictions  more  than 
one  day  in  advance  with  the  regression  model,  we  would  be  required  to  model  rainfall,  a difficult 
task.  Also,  as  Figures  26-30  indicate,  rainfall  has  different  affects  on  power  usage  at  different  times 
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of  the  year.  However,  we  are  not  able  to  easily  incorporate  this  observation  into  a regression  model 
framework.  For  these  reasons,  we  decided  to  develop  a base  model  which  could  be  used  as  a starting 
point  for  a new  year  and  whose  coefficients  could  be  updated  as  more  information  became  available. 
Once  the  base  model  had  been  decided  upon,  we  would  then  need  to  model  its  parameters.  Having 
accomplished  this  task,  we  would  be  able  to  make  predictions  for  the  days  remaining  based  on  the 
information  available  at  any  previous  point  in  time. 

Our  first  step  in  developing  a base  model  was  to  form  a theoretical  model  which  was  suitable 
for  all  of  the  available  years  (1992-1996).  We  returned  to  modeling  the  daily  maximum  series  {Yt} 
instead  of  the  relative  differenced  series  {X*},  because  the  peaks  of  {X*}  did  not  always  correspond 
with  the  peaks  of  {Y*}  and  our  main  concern  was  to  accurately  model  the  peaks  of  the  series  of 
daily  maxima.  We  began  our  search  for  a base  model  by  examining  sinusoidal  models  because  the 
series  {Yt}  resembled  a combination  of  sines  and  cosines. 

4.1  Sinusoidal  Models 

Only  the  first  three  available  years  (1992-1994)  were  considered  with  somewhat  discouraging  results. 
We  fit  the  sinusoidal  models  using  absolute  deviation  regression  (L\  regression).  This  procedure  fits 
a regression  model  by  minimizing  the  sum  of  the  absolute  deviations  of  the  errors  instead  of  the  sum 
of  the  squared  errors.  This  type  of  loss  function  was  chosen  in  an  effort  to  obtain  predictions  that 
were  more  robust  than  those  found  using  a squared  error  loss  function  ( L2  regression).  Figures  31- 
33  depict  both  the  observed  (‘o’)  and  predicted  (‘*’)  values  for  the  fitted  models  given  below. 


1992: 


7 rf  7 rf  jrf  7 rf 

Yt  = 3724.98-  788.11  cos(-)  - 1288.28cos(— ) - 1641.94cos(— ) + 428.04sin(— ) 
+383.80  sin(^)  - 2024.10  sinA  + 1999.57  sin(+  - 3422.79  sin(^) 


1993: 


7T  t 7r  f jrt  TTt 

Y,  = 2931.96  — 57.70  cos(—)  — 1201.66  cos(—)  - 949.12  cos(—)  + 88.89  sin(—) 

8 30  85  15 

+790.69 sin(— ) — 681.95  sin(^)  - 2203.08 sin(^) 


1994: 

TTt  7r  f 7rt  TCt 

Yt  = 4140.73  - 151.29  cos(-)  — 878.33  cos(-)  + 593.62 cos(— ) - 2733.19 cos(——) 

19  30  50  105 

+178.98 sin(j)  - 2209.81sin(^)  + 708.34sin(^) 
t = 1,2, ...,183. 

Only  two  of  the  predictor  variables  were  the  same  for  all  years,  and  minor  changes  in  the 
explanatory  variables  greatly  varied  the  fit.  This  made  it  difficult  to  arrive  at  a base  model  that 
would  be  suitable  for  all  of  the  years.  Therefore,  this  type  of  model  received  no  further  consideration. 
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o.'Y(t)  *:Yhat(l)  o:Y(t)  *:Yhat(t) 

2000  4000  6000  8000  10000  12000  0 2000  4000  6000  8000  10000  12000 


Figure  31:  Sinusoidal  Model  - 1992 


Figure  32:  Sinusoidal  Model  - 1993 
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Figure  33:  Sinusoidal  Model  - 1994 


4.2  Polynomial  Models 

We  next  considered  the  use  of  polynomial  functions  as  our  predictors  in  the  base  model.  The 
polynomials  we  decided  upon  were  shifted  Legendre  polynomials. 

4.2.1  Shifted  Legendre  Polynomials 

One  of  the  reasons  these  types  of  polynomials  were  considered  is  because  they  are  orthogonal  on 
the  interval  from  0 to  1.  The  first  ten  shifted  Legendre  polynomials  are  given  below. 

/o(z)  = 1 

fi  (x)  = 2x  — 1 

/2(x)  = 6x2  - 6x  + 1 

/3(x)  = 20x3  - 30x2  + 12x  - 1 

fi  (x)  = 70x4  - 140x3  + 90x2  - 20x  + 1 

h (x)  = 252x5  - 630x4  + 560x3  - 210x2  + 30x  - 1 

/is (x)  = 924x6  - 2, 772x5  + 3, 150x4  - 1, 680x3  + 420x2  - 42x  + 1 

fi (x)  = 3, 432x7  - 12, 012x6  + 16, 632x5  - 11, 550x4  + 4, 200x3  - 756x2  + 

56x  - 1 

/8(x)  = 12,  870x8  - 51, 480x7  + 84, 084x6  - 72,  072x5  + 34,  650x4  - 
9,  240x3  + 1, 260x2  - 72x  + 1 
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f9(x)  = 48,  620z9  - 218, 790a:8 + 411,  840a:7  - 420,420a:6  + 252, 252z5  - 
90,  090a:4  + 18,  480a:3  - 1,  980a:2  + 90a:  - 1 

Although  we  originally  used  these  ten  polynomials,  we  later  scaled  them  in  order  to  obtain  an 
orthonormal  set  of  polynomials  {p{}  (fj  pl(x)pj{x)dx  = 0 for  i ± j;  f*p?(x)dx  = 1).  The  scaling 
factors  used  to  ensure  orthonormality  for  the  ten  shifted  Legendre  polynomials  given  above  were 
1,  \/3,  \/5,  %/7 , 3,  \/TT,  vT3,  \/l5,  y/l7,  \/T9,  respectively.  The  orthonormal  property  is  nice 
because  it  simplifies  many  of  the  calculations  involving  the  design  matrix.  The  estimates  given  in 
this  section  are  based  on  the  scaled  polynomials  {p,- } . 

Models  consisting  of  various  combinations  of  the  scaled  polynomials  with  x = of  }^ear-90 
were  initially  fit  using  both  absolute  deviation  regression  (Li  regression)  and  least-squares  regression 
(Z/2  regression).  The  predictions  of  the  peaks  were  more  accurate  using  L\  regression  for  all  of  the 
years.  This  was  to  be  expected  since  the  severe  drops  in  usage  have  less  of  an  impact  under  this 
loss  function. 

The  model  which  was  the  most  parsimonious  without  losing  predictive  power  is  given  below. 

Y M = fioPo  {t)  + PiPi  ( t ) + P2P2  (0  + P3P3  {t)  + PaPs  (0  + Pz>P7  ( t ) + 

PeP9{t)  + e{t) 

• where 

t _ _1_  183 

183’'"’  183‘ 

Table  4.2.1  contains  the  parameter  estimates  and  standard  deviations  under  this  model  for  the  L\ 
regression  fit  for  the  years  1992-1996.  Figures  34-38  depict  both  the  observed  (‘o’)  and  predicted 
(‘*’)  values  for  these  years  using  this  model. 

The  use  of  shifted  Legendre  polynomials  resulted  in  a more  parsimonious  model  and  easier 
identification  of  each  polynomial’s  contribution  to  the  overall  fit  than  did  some  of  the  other  types 
of  orthogonal  polynomials  which  were  considered.  The  advantage  of  identifying  each  polynomial’s 
contribution  will  be  discussed  in  more  detail  in  Section  5.1. 

Table  4.2.1:  Parameter  Estimates  Using  Shifted  Legendre 
Polynomials:  L\  Regression 


Parameter 

1992 

1993 

1994 

1995 

1996 

Po 

3437.28 

(51.34) 

2820.11 

(48.41) 

4267.23 

(61.12) 

3946.93 

(70.38) 

3145.48 
( 66.98) 

Pi 

861.62 
( 51.31) 

1045.87 

(48.39) 

1219.31 

(61.09) 

1834.97 

(70.34) 

679.90 
( 66.94) 

Pi 

-1295.31 
( 51.50) 

-1067.69 

(48.56) 

-2295.49 

(61.31) 

-1089.79 

(70.60) 

-1098.90 

(67.19) 

P3 

-940.87 

(51.31) 

-1123.96 

(48.38) 

-1175.28 
( 61.09) 

-2097.92 

(70.34) 

-1290.70 
( 66.94) 

Pa 

957.84 

(51.26) 

617.87 

(48.34) 

576.30 

(61.03) 

529.66 

(70.28) 

890.92 

(66.88) 

p5 

-1033.36 

(51.16) 

-350.58 

(48.24) 

377.70 

(60.90) 

-300.39 

(70.13) 

-362.90 

(66.74) 

Pe 

563.55 

(50.96) 

344.79 

(48.05) 

-29.59 

(60.67) 

464.82 

(69.86) 

224.52 

(66.48) 
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Figure  34:  L\  Legendre  Polynomial  Regression  - 1992 


Figure  35:  L\  Legendre  Polynomial  Regression  - 1993 
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o:  Y(t)  *:  Yhat(t)  o:  Y(t)  *:  Yhat(t) 

2000  4000  6000  8000  10000  12000  0 2000  4000  6000  8000  10000  12000 


Figure  36:  L\  Legendre  Polynomial  Regression  - 1994 


Figure  37:  L\  Legendre  Polynomial  Regression  - 1995 
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Figure  38:  L\  Legendre  Polynomial  Regression  - 1996 


5 Improving  the  Base  Model 

Having  decided  which  orthogonal  polynomials  to  use  in  the  base  model,  we  began  to  refine  the 
model.  We  started  by  trying  to  obtain  more  accurate  predictions  of  the  peaks  in  the  data  by 
considering  different  types  of  loss  functions. 

5.1  Asymmetric  Loss  Functions 

Although  predictions  of  the  peaks  were  better  using  L\  regression  than  they  were  using  L2  regres- 
sion, they  were  still  not  satisfactory.  The  problem  with  using  this  type  of  loss  function  was  that  all 
of  the  extreme  points  received  less  weight,  including  the  peaks.  Therefore,  we  began  to  explore  the 
use  of  asymmetric  loss  functions,  functions  that  give  more  weight  to  positive  residuals.  Use  of  this 
type  of  loss  function  results  in  more  accurate  predictions  of  the  peaks  (as  well  as  poorer  predictions 
of  the  valleys,  but  this  was  not  of  concern  to  us). 

After  many  trials,  the  loss  function  we  decided  to  use  at  this  stage  was  an  absolute  deviation 
loss  which  gave  double  weight  to  the  positive  residuals.  Using  this  loss  function,  we  fit  the  same 
trend  model  as  in  Section  4.2.1.  This  model  is  reproduced  below. 

y M = PoPott)  + PlPl{t)  + p2P2(t)  + PsP3(t)  + teM  + P&7(t)  + 

P6P9{t)  + e(t)} 

where: 

Pi(t)  = the  ith  scaled,  shifted  Legendre  polynomial  at  time  t,  and 
t _ _1_  183 

183’ * ’ ’ ’ 183’ 
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The  parameter  estimates  for  this  model  are  given  in  Table  5.1.  Along  with  the  estimates  for  each 
year,  the  average  of  all  five  years  is  included  in  the  table.  With  the  exception  of  1994,  the  estimates 
for  each  year  were  of  the  same  sign  and  close  to  the  same  magnitude.  This  was  the  property  we 
had  been  searching  for  in  creating  a base  model. 


Table  5.1:  Parameter  Estimates  Based  on  an  Asymmetric  Loss  Function 


Parameter 

1992 

1993 

1994 

1995 

1996 

average 

Po 

3855.73 

3385.58 

4682.06 

4423.81 

3833.92 

4036.22 

Pi 

1069.31 

1302.74 

1317.00 

1878.50 

820.47 

1277.61 

P2 

-1378.47 

-1385.32 

-2204.87 

-1020.76 

-1449.68 

-1487.82 

Ps 

-1125.87 

-1383.07 

-1216.69 

-2335.39 

-1561.61 

-1524.53 

04 

1018.92 

747.14 

425.41 

507.56 

1145.80 

768.97 

(k 

-1079.67 

-376.33 

332.89 

-427.71 

-492.58 

-408.68 

Pe 

509.18 

358.21 

-30.30 

455.44 

293.01 

317.11 

Plots  of  the  observed  values,  ‘o’,  and  both  sets  of  the  predicted  values  (those  using  the  parameter 
estimates  for  the  year,  and  those  using  the  average  parameter  estimates,  ‘A’)  are  given  in 
Figures  39  through  43.  Use  of  this  loss  function  resulted  in  more  accurate  predictions  of  the  peaks 
than  did  the  previously  considered  symmetric  loss  functions.  Additionally,  with  the  understanding 
that  the  estimates  were  not  meant  to  be  used  as  predictions,  simply  as  our  best  guess  for  a new 
year,  the  average  estimates  were  a promising  start  for  our  base  model. 

Plots  of  the  polynomials  p,-,  i=l,2,3,5,7,9,  times  their  estimated  coefficients  /3j,  j=l,...,6,  are 
given  in  Figures  44  through  49.  These  graphs  show  the  contribution  of  each  polynomial  to  the 
overall  fit  of  the  model.  With  a few  exceptions,  each  polynomial’s  contribution  is  similar  from  one 
year  to  the  next.  Again,  this  reinforced  our  belief  that  we  had  found  an  adequate  base  model. 
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Figure  39:  Asymmetric  Loss  Function-  1992 


Figure  40:  Asymmetric  Loss  Function-  1993 
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Figure  41:  Asymmetric  Loss  Function-  1994 


Figure  42:  Asymmetric  Loss  Function-  1995 
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o:  Y(t)  *:  Yhat(t)  A:  averageYhat(t) 


Figure  43:  Asymmetric  Loss  Function-  1996 


Figure  44:  (3\  * pi 
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Figure  45:  (32  * p2 


P3 

Figure  46:  /33  * p3 
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Figure  47:  (34*  p5 


Figure  48:  (3$  * p 7 
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Figure  49:  (36  * p9 


5.2  Integrating  Rainfall 

It  has  already  been  shown  that  rainfall  plays  an  important  role  in  predicting  maximum  daily 
power  usage.  One  of  the  ways  we  have  considered  including  this  information  is  to  use  rainfall  as  a 
predictor  for  the  residuals  obtained  under  the  polynomial  fit.  However,  this  method  only  improved 
the  predictions  of  the  drops  in  usage,  not  the  peaks.  Therefore,  the  benefit  of  using  rainfall  in  this 
manner  was  not  worth  the  required  effort  of  predicting  rainfall. 

We  also  considered  including  rainfall  as  one  of  the  predictors  in  the  base  model.  However,  the 
predictions  under  this  modeling  scheme  were  not  as  accurate,  the  model  was  more  complicated, 
and  the  range  of  magnitude  of  each  parameter  was  greater  from  one  year  to  the  next.  At  this 
point,  we  decided  to  postpone  the  inclusion  of  rainfall  and  move  on  to  the  next  step,  modeling  the 
parameters. 

6 Model  Choices  for  the  Parameters  of  the  Base  Model 

We  wanted  to  update  the  parameters  of  the  base  model  as  more  information  became  available.  The 
structural  model  framework  allows  the  parameters  to  change  with  time  and  makes  it  easy  to  obtain 
predictions  of  the  power  usage  for  the  remaining  days  based  on  any  previous  day.  Therefore,  we 
began  modeling  the  parameters  according  to  the  structural  model  given  below. 

Yt  = p,/3 , + Wtl  W(~N(0,<ri) 

0t+ 1 = 0t  + V(,  V(  ~ N(0,<3) 
t = 1,  — , 183, 
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where  : 


Yt 

P t 

a 

with  initial  value 
/3i 
and 


Q 


the  maximum  daily  power  usage  in  kw 
Po(t)  Pl(t)  p2{t)  P3(t)  Ps(t)  p7(t)  p9{t) 

Po(t)  Pi(t)  h(t)  p3{t)  Pa  (t)  Ps  {t)  Pe{t) 


4036  1278  -1488  -1525  769  -409  317 

al  0 0 0 0 0 0 

0 a\  0 0 0 0 0 

0 0 o\  0 0 0 0 

0 0 0 a\  0 0 0 

0 0 0 0 a\  0 0 

0 0 0 0 0 erg  0 

0 0 0 0 0 0 a\ 


Under  this  modeling  scheme,  the  Kalman  recursions  simplify  as  follows: 

Pt+i  = Pt  + At_1  {Yt  ~ P tPt) 

Qt+i  — fit  + Q — 0*Af  *0j 

where 

At  = pttttp't  + vl,  and 

0*  = tttp't- 

For  the  initial  conditions,  we  let  (3l  equal  the  average  of  the  parameter  estimates  obtained  in 
Section  5.1  and  The  best  linear  mean-squared  h-step  predictors  PtYt+h  are  determined 

by  the  following  equations: 

PtPt+h  = PtPt+i  =Pt+u  and 

PtYt+h  — Pt+Zi^t+i* 

With  this  model,  to  obtain  a prediction  for  day  s (for  example)  based  on  the  information 
available  at  time  t (s  > t),  the  interested  party  would  simply  calculate  ps(3t+l.  Therefore,  the 
prediction  for  the  remaining  days  based  on  the  information  available  at  time  t would  be  a polynomial 
fit  similar  to  the  base  model  but  with  different  coefficients. 

The  above  model  was  fit  using  a Fortran  program  similar  to  the  one  used  in  Section  3.1.  The 
function  that  the  program  seeks  to  optimize  has  many  local  minimums.  Therefore,  these  estimates 
are  not  obtained  from  a function  that  has  been  minimized  globally,  hence  they  are  not  uniqe 
and  may  not  be  the  true  maximum  likelihood  estimates.  However,  reasonable  variations  of  these 
estimates  result  in  almost  identical  predictions.  Thus,  we  did  not  feel  it  was  necessary  to  spend  the 
time  required  to  find  the  exact  maximum  likelihood  estimates  if,  in  fact,  they  exist. 

In  our  initial  fit,  we  assumed  that  Q = ctqI7.7  (<Jq  = a\  = . . . = eg).  The  maximum  likelihood 
estimates  of  the  parameters  under  this  model  are  given  in  Table  6.  Figures  50  through  54  depict 
the  observed  values  (‘o’)  and  the  one-step  predictions  (‘*’)  for  this  model.  Comparison  of  these 
figures  with  those  resulting  from  earlier  models  shows  the  gain  in  predictive  power  achieved  using 
the  base  model  approach. 
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Next  we  fit  the  same  theoretical  model  but  allowed  the  variances  of  the  /Ps  to  differ  (o-q  ^ o\ 

. . . ^ (Jq).  Removal  of  the  equal  variances  constraint  resulted  in  aj  through  <r|  being  close  to  zero 
and  <7q  being  quite  large.  Therefore,  the  one-step  predictions  were  essentially  found  by  moving  the 
base  model  polynomial  up  or  down  so  that  it  passed  through  the  observed  data.  This  resulted  in 
very  poor  h- step  predictions  (h  > 1).  For  this  reason,  we  kept  the  constraint  of  equal  variances  in 
our  model  framework.  The  interest  at  this  point  was  in  improving  the  predictions  obtained  under 
this  model. 

Table  6:  Maximum  Likelihood  Estimates 


Parameter 

1992 

1993 

1994 

1995 

1996 

29,108 

86,140 

45,242 

41,222 

230,771 

136,514 

142,993 

125,771 

77,909 

212,309 

Figure  50:  One-step  Predictions  from  Structural  Model  - 1992 
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Figure  51:  One-step  Predictions  from  Structural  Model  - 1993 


Figure  52:  One-step  Predictions  from  Structural  Model  - 1994 
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Figure  53:  One-step  Predictions  from  Structural  Model  - 1995 


Figure  54:  One-step  Predictions  from  Structural  Model  - 1996 
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7 Improving  the  Parameter  Modeling  Strategy 

7.1  Robust  Kalman  Filter 


Although  the  parameters  of  the  base  model  were  fit  using  a robust  technique,  the  “updated” 
parameter  estimates,  (3t,  from  Section  6 obtained  using  the  Kalman  recursions  were  not  robust. 
This  is  because  the  Kalman  equations  give  the  best  linear  mean-square  predictors  which  are  heavily 
influenced  by  outlying  observations.  At  this  point  we  began  researching  alternative  ways  to  arrive 
at  robust  parameter  estimates. 

Our  search  resulted  in  the  discovery  of  a paper  entitled,  Robust  Kalman  Filter  and  Its  Applica- 
tion in  Time  Series  Analysis,  by  Tomas  Cipra  and  Rosario  Romera  (Cipra  and  Romera  [6]).  This 
paper  suggested  converting  the  structural  model  to  a linear  regression  model  and  then  minimizing 
a suitable  robustifying  function  in  order  to  obtain  a robust  form  of  the  Kalman  recursions. 

Using  these  results,  we  began  searching  for  a suitable  robustifying  function.  We  first  tried  what 
is  known  as  the  Huber  function.  This  function  is  given  below. 


h{z ) = 


for  \z\  < c, 
csgn(z),  for  \z\  > c, 


where  c is  some  constant.  (The  reader  is  referred  to  Lehmann  [9]  or  Huber  [8]  for  further  information 
on  the  Huber  function.)  This  function  did  not  provide  us  with  satisfactory  predictions  because  it 
gives  equal  weight  to  positive  and  negative  residuals  that  are  of  the  same  magnitude. 

We  eventually  decided  upon  an  asymmetric  squared  error  function  which  weights  the  positive 
residuals  ten  times  greater  than  the  negative  ones.  The  number  ten  is  somewhat  arbitrary.  As 
long  as  the  positive  residuals  receive  significantly  more  weight,  the  fit  will  be  essentially  the  same. 
Applying  the  results  of  the  paper  by  Cipra  and  Romera  to  this  function,  the  robust  recursive 
formulas  for  obtaining  /3t+1  are  as  follows: 


Pt+i 

n*+i 

where 


Pt  + 


Fit  ~\-  Q — 


pAp!  + oi /c 

0*©; 


(Yt  - P A) 


PtfitPt  + ffi/c’ 


10,  if  {Yt  - Yt)  > 0, 
1,  otherwise. 


For  the  case  when  Yt  < Yt,  we  obtain  the  same  predictors  as  in  Section  6. 

In  order  to  simplify  the  optimization  process,  we  chose  <7^  = 10,  000  and  Oq  = 100  (recall  Q = 
Cq  1-1.7).  Again,  these  numbers  were  chosen  somewhat  arbitrarily.  Although  several  optimization 
programs  were  used  to  find  these  estimates,  the  function  has  many  local  minimums  making  it 
difficult  to  arrive  at  a global  minimum.  However,  the  fit  changed  very  little  for  reasonable  variations 
of  these  estimates.  This  was  beneficial  not  only  for  the  reason  just  mentioned,  but  also  because  we 
were  able  to  use  the  same  estimates  for  all  five  years. 


7.2  Incorporating  Weather  Variables 

We  began  including  weather  variables  in  our  modeling  framework  at  the  same  time  as  our  search  for 
more  robust  techniques  was  initiated.  Rainfall  was  the  first  weather  variable  considered.  We  knew 
that  rainfall  was  responsible  for  the  extreme  drops  in  usage  and  that  large  amounts  of  precipitation 
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had  a prolonged  affect  in  many  cases.  We  wanted  to  use  this  information  to  prevent  the  parameter 
estimates  from  changing  drastically  in  these  cases,  not  to  obtain  more  accurate  predictions  of  the 
drops.  Therefore,  we  allowed  the  variance  of  the  polynomial  coefficients  ( ctq ) to  change  according 
to  rainfall.  This  was  done  by  first  creating  a function  to  account  for  the  additive  affect  of  rainfall. 
After  much  debate  and  experimentation,  we  settled  on  the  function  h(t)  = 0.6 R(t)  + QAh(t  - 1), 
where  R(t)  was  the  rainfall  for  day  t.  We  let  a%(t)  = 100exp-3O14/l(*).  This  technique  allows  us  to 
include  the  effects  of  rainfall  in  our  one-step  predictions  without  having  to  predict  precipitation. 
For  our  h- step  predictions  (h  > 1),  we  predicted  rainfall  to  be  0mm  with  satisfying  results.  The 
value  3.014  was  found  by  averaging  the  optimum  values  of  a in  the  function  100  exp“a/lW  for  all 
five  available  years.  The  optimum  values  of  a were  found  using  the  IMSL  routine  DUVMIF  and 
are  given  in  Table  7.2. 


Table  7.2:  Optimum  values  for  a 


1992 

1993 

1994 

1995 

1996 

average 

a 

4.475 

-0.567 

2.049 

1.719 

7.404 

3.014 

In  addition  to  including  rainfall,  we  investigated  the  inclusion  of  ETAcorn.  We  began  by  adding 
it  to  our  model  in  the  same  manner  as  we  had  for  precipitation.  However,  this  did  not  significantly 
improve  our  predictions.  We  then  included  ETAcorn  as  a parameter  in  the  base  model.  Although 
this  requires  the  user  to  have  an  estimate  of  ETAcorn  for  the  entire  year,  it  can  be  calculated  at  the 
beginning  of  the  year  assuming  normal  conditions.  Including  ETAcorn  in  this  manner  did  improve 
our  predictions.  The  resulting  model  (as  given  below)  is  a variation  of  that  given  in  Section  6. 


Yt  = 

Pt+ i — 
t = 
where  : 

Yt  = 
P t = 


Pt  = 

with  initial  value 

Pi  = 
°l(t)  = 

and 

hit)  = 


PtA  + Wi,  Wt  ~ N(0, 10,  000) 
/3t  + Vt,  V<  ~ N(0,  0q(£)  • h-&) 
!,•••,  183, 


the  maximum 

daily  power  usage  in 

kw 

Po(t) 

Pl(t ) 

P2<«) 

Ps(t) 

Ps(t) 

PtW 

Ps(t) 

eta[t)  j 

’ 0o(t) 

m 

Pl(t) 

AW 

AW 

AW 

AW 

/M0  ] 

3793  1083  -1167  -1344  637  -318  245  268 
100exp-3-ol4,/l(t) 

0.6R{t)  + 0.4h(t-  1). 


For  the  initial  conditions,  we  let  equal  an  8 by  8 matrix  of  zeros  and  we  chose  fli  to  be  the 
average  of  the  parameter  estimates  obtained  in  a similar  manner  to  those  found  in  Section  5.1.  The 
asymmetric  loss  function  used  in  obtaining  /3\  for  this  model  was  ten  times  the  squared-error  loss  for 
positive  residuals  and  squared-error  loss  for  negative  residuals.  The  h- step  predictions  are  found  in 
the  same  manner  as  that  described  in  Section  6.  Namely,  the  prediction  of  maximum  power  usage 
for  day  s based  on  the  information  available  up  to  time  t {t  < s)  is  found  by  calculating  ps(3t+i . 
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Five  sets  of  plots  depicting  the  observed  values  (‘o’)  and  h-step  predictions  (**’),  h = 1,3, 5,  7, 
for  all  years  under  this  model  are  included  below  as  Figures  55-74  ( Figures  55-59  give  the  one-step 
predictors,  Figures  60-64  the  three-step  predictors,  Figures  65-69  the  five-step  predictors  and  Fig- 
ures 70-74  the  seven-step  predictors).  The  plotted  predicted  values  (**’)  are  actually  the  maximum 
of  zero  and  the  values  predicted  by  the  model.  The  plots  of  the  one-step  predictors  include  darkened 
circles  indicating  that  there  was  a positive  amount  of  precipitation  on  the  previous  day.  As  was 
discussed  earlier,  the  polynomial  coefficients  will  change  by  a negligible  amount  when  there  has 
been  a significant  amount  of  rainfall;  thus  the  predictions  will  closely  follow  a polynomial  curve  for 
the  next  several  days  after  a large  rainfall.  This  effect  due  to  rainfall  can  be  seen  by  careful  exam- 
ination of  Figures  55-59.  Comparison  of  the  one-step  predictions  from  this  model  with  predictions 
from  earlier  models  shows  the  gain  in  predictive  power  obtained  by  use  of  this  modeling  scheme. 
Even  the  predictions  made  seven  days  in  advance  are  reasonably  good  when  one  considers  that 
they  are  based  only  on  historical  data,  past  rainfall  and  ETAcorn.  Residual  plots  for  the  one-step 
and  seven-step  predictors  are  given  in  Figures  75-84.  Residuals  are  calculated  as  the  observed  value 
minus  the  predicted  value. 

The  Fortran  program  used  to  fit  the  model  described  in  this  section  is  similar  to  the  programs 
used  in  Sections  3.1  and  6 with  adjustments  made  to  incorporate  the  changes  discussed  above.  In 
addition,  the  use  of  Hooke  and  Jeeve’s  algorithm  to  find  the  maximum  likelihood  estimates  of  cr^ 
and  Q is  no  longer  required  since,  for  this  model,  crj  is  fixed  and  Oq  is  determined  by  the  function 
h(t)  previously  defined. 


43 


o:Y(t)  *:  1-step  o:  Y(t)  *:  1-step 

)0  6000  8000  10000  12000  0 2000  4000  6000  8000  10000  12000 


Figure  55:  One-step  Predictions  from  Final  Model  - 1992 


Figure  56:  One-step  Predictions  from  Final  Model  - 1993 
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o:  Y(t)  *:  1-step  o:  Y(t)  *:  1-step 


Figure  57:  One-step  Predictions  from  Final  Model  - 1994 


Figure  58:  One-step  Predictions  from  Final  Model  - 1995 
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o:  Y(t)  *:  3-step  o:  Y(t)  *:  1-step 


Figure  59:  One-step  Predictions  from  Final  Model  - 1996 


Figure  60:  Three-step  Predictions  from  Final  Model  - 1992 
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o:  Y(t)  *:  3-step  o:  Y(t)  *:  3-step 


Figure  61:  Three-step  Predictions  from  Final  Model  - 1993 


Figure  62:  Three-step  Predictions  from  Final  Model  - 1994 
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o:  Y(t)  *:  3-step  o:  Y(t)  *:  3-step 


Figure  63:  Three-step  Predictions  from  Final  Model  - 1995 


Figure  64:  Three-step  Predictions  from  Final  Model  - 1996 
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o:  Y(t)  *:  5-step  o:Y(t)  *:  5-step 


Figure  65:  Five-step  Predictions  from  Final  Model  - 1992 


Figure  66:  Five-step  Predictions  from  Final  Model  - 1993 
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o:  Y(t)  *:  5-step  o:Y(t)  *:  5-step 


Figure  67:  Five-step  Predictions  from  Final  Model  - 1994 


Figure  68:  Five-step  Predictions  from  Final  Model  - 1995 
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o:  Y(t)  *;  7-step  o:Y(t)  *:  5-step 
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Figure  69:  Five-step  Predictions  from  Final  Model  - 1996 


Figure  70:  Seven-step  Predictions  from  Final  Model  - 1992 
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o:  Y(t)  *:  7-step  o:Y(t)  *:  7-step 


Figure  71:  Seven-step  Predictions  from  Final  Model  - 1993 


Figure  72:  Seven-step  Predictions  from  Final  Model  - 1994 
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o:  Y(t)  *:  7-step  o:  Y(t)  *:  7-step 


Figure  73:  Seven-step  Predictions  from  Final  Model  - 1995 


Figure  74:  Seven-step  Predictions  from  Final  Model  - 1996 
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Figure  75:  Residuals  from  One-step  Predictions  - 1992 


Figure  76:  Residuals  from  One-step  Predictions  - 1993 
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Figure  77:  Residuals  from  One-step  Predictions  - 1994 
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Figure  78:  Residuals  from  One-step  Predictions  - 1995 
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Figure  79:  Residuals  from  One-step  Predictions  - 1996 
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Figure  80:  Residuals  from  Seven-step  Predictions  - 1992 
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Figure  81:  Residuals  from  Seven-step  Predictions  - 1993 
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Figure  82:  Residuals  from  Seven-step  Predictions  - 1994 
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Figure  83:  Residuals  from  Seven-step  Predictions  - 1995 
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Figure  84:  Residuals  from  Seven-step  Predictions  - 1996 
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8 Modeling  New  Data 

Although  the  focus  of  this  project  was  on  the  power  data  obtained  from  the  Holyoke  station,  our 
goal  was  to  develop  a model  which  could  be  used  to  make  power  predictions  for  other  stations  in 
the  surrounding  area.  In  this  section  we  examine  the  predictions  under  the  final  model  (described 
in  Section  7.2)  for  two  additional  power  stations,  Haxtun  and  Interstate. 

Figures  85-89  give  the  one-step  predictors  of  the  maximum  daily  power  usage  for  the  Haxtun 
station,  Figures  90-94  give  the  one-step  predictors  for  the  Interstate  station,  Figures  95-99  give  the 
seven-step  predictors  for  the  Haxtun  station  and  Figures  100-104  give  the  seven-step  predictors  for 
the  Interstate  station.  In  addition,  each  graph  depicts  the  observed  daily  maxima  for  the  station 
under  consideration.  The  time  frame  used  was  the  same  as  that  for  the  Holyoke  station,  namely 
April  1 to  September  30.  In  order  to  show  the  effects  of  rainfall,  solid  circles  have  been  plotted 
indicating  a positive  amount  of  precipitation  on  the  previous  day.  Although  the  magnitude  of  the 
daily  maxima  was  not  the  same  for  these  stations  (the  maxima  for  the  Haxtun  station  are  somewhat 
smaller  than  for  Holyoke  whereas  the  maxima  for  the  Interstate  station  are  roughly  double),  the 
model  adapted  quickly  and  the  predictions  were  quite  good.  For  the  first  three  years  it  appeals 
as  if  missing  data  at  the  beginning  of  the  year  was  filled  in  with  zeros  for  the  Interstate  station 
(see  Figures  90-92).  Even  in  this  situation,  the  predictions  for  the  first  several  days  of  true  data 
were  satisfactory,  and  the  ability  of  the  model  to  adjust  to  the  sudden  change  (due  both  to  the 
missing  observations  and  the  magnitude  of  the  observations)  is  amazing.  Within  just  a few  days, 
the  estimates  of  the  polynomial  coefficients  adapted  to  adjust  for  this  change. 


Figure  85:  One-step  Predictors  for  the  Haxtun  Station  - 1992 
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o:  Y(t)  *:  1-step  o:  Y(t)  *:  1-step 
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Figure  86:  One-step  Predictors  for  the  Haxtun  Station  - 1993 


Figure  87:  One-step  Predictors  for  the  Haxtun  Station  - 1994 


60 


Figure  88:  One-step  Predictors  for  the  Haxtun  Station  - 1995 


Figure  89:  One-step  Predictors  for  the  Haxtun  Station  - 1996 
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o:Y(t)  *:  1-step  o:  Y(t)  *:  1-step 

5000  10000  15000  20000  0 5000  10000  15000  20000 


Figure  90:  One-step  Predictors  for  the  Interstate  Station  - 1992 


Figure  91:  One-step  Predictors  for  the  Interstate  Station  - 1993 
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o:Y(t)  *:1-step  o:  Y(t)  *:  1-step 
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Figure  92:  One-step  Predictors  for  the  Interstate  Station  - 1994 
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Figure  93:  One-step  Predictors  for  the  Interstate  Station  - 1995 
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o:  Y(t)  *:  7-step  o:  Y(t)  *:  1-step 
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Figure  94:  One-step  Predictors  for  the  Interstate  Station  - 1996 


Figure  95:  Seven-step  Predictors  for  the  Haxtun  Station  - 1992 
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o:  Y(t)  *:  7-step  o:  Y(t)  *:  7-step 


Figure  96:  Seven-step  Predictors  for  the  Haxtun  Station  - 1993 


Figure  97:  Seven-step  Predictors  for  the  Haxtun  Station  - 1994 
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o:  Y(t)  *:  7-step  o:  Y(t)  *:  7-step 


Figure  98:  Seven-step  Predictors  for  the  Haxtun  Station  - 1995 


Figure  99:  Seven-step  Predictors  for  the  Haxtun  Station  - 1996 
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o:  Y(t)  *:  7-step  o:Y(t)  *:  7-step 
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Figure  100:  Seven-step  Predictors  for  the  Interstate  Station  - 1992 


Figure  101:  Seven-step  Predictors  for  the  Interstate  Station  - 1993 
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o:  Y(t)  *:  7-step  o:  Y(t)  *:  7-step 
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Figure  102:  Seven-step  Predictors  for  the  Interstate  Station  - 1994 


Figure  103:  Seven-step  Predictors  for  the  Interstate  Station  - 1995 
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Figure  104:  Seven-step  Predictors  for  the  Interstate  Station  - 1996 

In  order  to  best  utilize  the  current  model,  we  suggest  two  prediction  strategies.  The  first, 
Strategy  A,  is  to  be  used  in  the  situation  when  predictions  must  be  made  at  the  beginning  of  the 
week  for  the  next  seven  days.  In  order  to  protect  against  under-predicting  the  peaks,  we  suggest 
that  the  predictions  for  the  next  seven  days  based  on  the  current  day  (e.g.  the  predictions  for  days 
i -f  1,  t + 2,  ...,£  + 7 based  on  day  t ) all  be  the  same  value.  This  value  is  to  be  the  maximum  of  the 
predictions  for  days  t + 1, . . . , t + 7 based  on  day  t,  day  t - 1 and  day  t — 2,  as  well  as  the  observed 
maximums  for  days  £,£  — 1 and  t - 2.  If  it  is  possible  to  run  the  program  daily,  we  recommend 
Strategy  B.  During  the  first  month  of  the  growing  season,  April,  Strategy  B takes  the  maximum  of 
the  predictions  for  days  t+ 1, . . . , t+ 7 based  on  day  t as  its  prediction  for  the  following  day.  During 
the  remaining  months,  May  thru  September,  the  prediction  for  day  t + 1 is  found  by  taking  the 
maximum  of  the  predictions  for  days  t -f-  1,  • • • , t -f  4 based  on  day  t and  the  observed  maximums 
for  days  £,  t — 1,  t — 2. 

Since  the  maximum  usage  each  month  was  of  particular  interest,  we  have  focused  on  the  ac- 
curacy of  its  prediction  in  the  table  and  graphs  that  follow.  Table  8 gives  the  predictions  of  each 
monthly  maximum  under  both  modeling  strategies,  A and  B,  for  the  Holyoke  data  for  all  five  avail- 
able years.  Also  included  are  the  1 and  7-step  predictors  discussed  previously.  Comparison  of  the 
methods  shows  that  while  Strategies  A and  B are  both  more  accurate  than  the  one  and  seven-step 
predictors,  Strategy  B is  to  be  favored.  This  is  again  apparent  in  comparison  of  Figures  105-107 
to  Figures  108-110.  The  first  set  of  bar  graphs  (Figures  105-107)  show  the  observed  and  predicted 
monthly  maximums  for  all  five  years  for  the  Holyoke,  Haxtun  and  Interstate  stations  under  Strategy 
A.  The  remaining  graphs  (Figures  108-110)  depict  the  same  under  Strategy  B.  Examination  of  the 
figures  shows  that  while  the  predictions  under  both  strategies  are  quite  good,  there  is  a significant 
amount  of  improvement  in  using  Strategy  B over  Strategy  A.  Considerations  on  how  to  further 
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improve  this  modeling  strategy  are  given  in  the  following  section. 


Table  8:  Monthly  Maximum  Predictions  for  Holyoke 


Year 

Month 

Observed 

1-step 

7-step 

Strategy  A 

Strategy  B 

1992 

April 

3070 

2408 

1989 

2436 

2676 

May 

4516 

3717 

2407 

2436 

3997 

June 

3531 

2871 

3089 

4529 

3717 

July 

9366 

9109 

8946 

9612 

9404 

August 

8450 

6178 

4874 

7603 

7972 

September 

4706 

4176 

3427 

4874 

5164 

1993 

April 

1146 

1240 

* 

1959 

1959 

May 

2189 

2002 

1841 

2043 

2178 

June 

7286 

6349 

3956 

4116 

7339 

July 

9084 

7741 

7741 

5219 

8778 

August 

8888 

8671 

5630 

5219 

9084 

September 

4539 

3627 

2472 

4429 

3911 

1994 

April 

1763 

1779 

1428 

1505 

2198 

May 

2926 

2475 

2696 

2788 

2559 

June 

9406 

9140 

10,021 

10,284 

9350 

July 

9538 

9069 

8430 

10,162 

10,674 

August 

8692 

7882 

7826 

8300 

8316 

September 

4470 

5375 

6097 

8052 

6296 

1995 

April 

1486 

1611 

* 

1959 

1709 

May 

1613 

1399 

1118 

1222 

1499 

June 

2310 

2255 

3590 

3933 

2984 

July 

9700 

8842 

5284 

6282 

8922 

August 

9308 

9178 

9866 

10,107 

9274 

September 

7770 

7512 

7896 

8842 

8398 

1996 

April 

3036 

2780 

2135 

2139 

2836 

May 

4320 

3286 

2572 

2663 

3629 

June 

5414 

4977 

3750 

4382 

5533 

July 

9740 

9412 

10,746 

10,719 

9534 

August 

8876 

7826 

6897 

7585 

8518 

September 

3012 

1064 

3096 

4514 

2517 

* the  monthly  maximum  occured  before  the  7th  day 
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Figure  105:  Monthly  Maximums  for  Holyoke  Under  Strategy  A 
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Figure  106:  Monthly  Maximums  for  Haxtun  Under  Strategy  A 
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Figure  107:  Monthly  Maximums  for  Interstate  Under  Strategy  A 
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Figure  108:  Monthly  Maximums  for  Holyoke  Under  Strategy  B 
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Figure  109:  Monthly  Maximums  for  Haxtun  Under  Strategy  B 
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Figure  110:  Monthly  Maximums  for  Interstate  Under  Strategy  B 
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9 Further  Considerations 

9.1  Wavelets 


One  of  the  directions  we  have  considered  taking  in  predicting  the  data  concerns  the  use  of  wavelet 
functions  in  lieu  of  Legendre  polynomials.  Wavelets  were  considered  because,  theoretically,  a small 
number  of  wavelets  can  be  used  to  closely  model  the  data  even  when  spikes  are  present.  Polynomial 
functions  tend  to  be  much  smoother  than  wavelets  which  makes  it  difficult  to  model  erratic  data 
in  detail  without  the  use  of  many  functions. 

A wavelets  module  is  available  for  S-Plus  which  we  used  in  our  preliminary  fitting  of  wavelet 
functions  to  the  data.  Applied  Wavelet  Analysis  with  S-PLUS  by  Bruce  and  Gao  ([3])  was  used 
els  a reference  to  this  module.  The  module  contains  several  functions  which  can  be  used  in  fitting 
wavelets  to  the  data.  We  tried  many  of  these  and  found  the  function  best.basis  to  be  the  most 
satisfactory.  Figures  111-115  depict  the  original  series  (dashed  line)  along  with  the  reconstructed 
series  (solid  line).  The  reconstruced  series  was  found  by  first  applying  the  best.basis  function  to  the 
original  series  using  the  dS  wavelet,  abstracting  the  top  ten  wavelet  functions  (those  containing  the 
most  information)  and  then  reconstructing  the  series  from  these  ten  wavelet  functions  (see  Bruce 
and  Gao  [3]  for  further  information).  However,  we  came  across  some  of  the  same  problems  we 
had  encountered  in  fitting  the  polynomial  model.  For  one,  the  robust  fit  discounted  all  extreme 
values,  large  and  small.  Also,  the  wavelet  functions  which  best  fit  the  data,  among  those  functions 
considered,  differed  from  one  year  to  the  next.  At  this  point,  we  decided  that  there  was  insufficient 
time  to  explore  these  problems  further.  However,  since  the  preliminary  plots  were  encouraging, 
wavelet  functions  are  mentioned  here  as  a possible  route  to  consider  in  further  studies  involving 
similar  data  sets. 


Figure  111:  Wavelet  Reconstruction  - 1992 
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Figure  112:  Wavelet  Reconstruction  - 1993 


Figure  113:  Wavelet  Reconstruction  - 1994 
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Figure  114:  Wavelet  Reconstruction  - 1995 


Figure  115:  Wavelet  Reconstruction  - 1996 
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9.2  Spatial  Analysis 

Another  direction  we  have  not  yet  explored  is  the  use  of  a spatial  factor  in  our  model.  At  this 
point,  only  the  weather  variables  from  one  of  the  “nearby”  stations  have  been  integrated  into  our 
modeling  scheme.  For  many  of  the  power  stations,  there  are  other  weather  stations  which  are 
closer  than  the  Holyoke  weather  station.  The  weather  data  collected  from  these  stations  could  be 
combined  to  provide  more  information  resulting  in  better  predictions.  Additionally,  the  distance 
between  the  power  station  and  the  Holyoke  weather  station  is  much  further  for  some  stations  than 
for  others.  This  could  be  taken  into  account  by  allowing  weather  variables  to  have  a diminishing 
affect  based  on  the  distance  between  the  weather  and  power  stations.  It  is  also  likely  that  the  daily 
maxima  for  power  stations  that  are  close  to  one  another  are  more  strongly  related  than  the  maxima 
for  stations  that  are  far  apart.  Again,  this  potential  relationship  could  be  included  in  our  model 
by  way  of  a spatial  component.  The  idea  of  a spatial  element  is  mentioned  here  as  a consideration 
for  further  analyses. 
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